/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of UHM (Unassembled HyperMatrix) sparse direct solver.
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/

#include "petscksp.h"

#include "uhm.hxx"
#include "uhm/interf/petsc.hxx"

#define UHM_ERROR_TOL 1.0e-5

int main (int argc, char **argv) {

  // input check
  if (argc != 5) {
    printf("Try : %s [n_thread][decomposition][is_spd][input_file]\n", argv[0]);
    return 0;
  }

  FLA_Init();

  uhm::Mesh m;
  uhm::Scheduler s;

  int n_threads, decomposition, is_spd;
  char *filename;

  n_threads     = atoi( (argv[1]) );
  decomposition = atoi( (argv[2]) );
  is_spd        = atoi( (argv[3]) );
  filename      = argv[4];

  uhm::set_num_threads(n_threads);
  uhm::set_hier_blocksize(256);

  printf( "BEGIN : Import mesh from file : %s\n", filename );
  m = new uhm::Mesh_;
  m->import_file( filename );
  printf( "END   : Import \n" );

  printf( "BEGIN : Build tree \n" );
  uhm::build_tree(m);
  printf( "END   : Build tree\n" );
  printf( "New mesh info < n_nodes %d, n_elements %d >\n", 
	  m->get_n_nodes(), m->get_n_elements() );
  
  m->lock();

  int datatype = UHM_REAL; 
  int n_rhs    = 1;

  printf( "BEGIN : Create matrix without buffer \n" );
  m->create_matrix_without_buffer( datatype, n_rhs );
  printf( "END   : Create matrix without buffer \n" );

  printf( "BEGIN : Create buffer \n" );
  m->create_leaf_matrix_buffer();
  printf( "END   : Create buffer \n" );

  switch (decomposition) {
  case UHM_CHOL: 
    m->random_spd_matrix(); 
    m->triangularize();
    break;
  default: 
    m->random_matrix(); 
    break;
  }

  m->set_rhs();

  // ----------------------------------------------------------------

  double 
    t_base, 
    t_init, 
    t_solve, 
    t_check,
    t_finalize;

  uhm::interf::Petsc petsc = new uhm::interf::Petsc_(UHM_REAL);

  // Step -1 : initialization
  t_base       = uhm::timer();
  petsc->init(argc, argv);
  t_init = uhm::timer() - t_base;

  petsc->set_show_n_rhs(10);

  // Step 0 : export matrix
  m->export_matrix(petsc, n_rhs, is_spd);

  petsc->create_ksp();
  //  petsc->set_pc_on_ksp(PCNONE);
  petsc->set_pc_on_ksp(PCBJACOBI);

  printf( "BEGIN : Solve \n" );
  t_base       = uhm::timer();
  petsc->solve();
  t_solve = uhm::timer() - t_base;
  printf( "END   : Solve \n" );

  t_base       = uhm::timer();
  petsc->check();
  t_check      = uhm::timer() - t_base;

  petsc->delete_ksp();
  petsc->delete_sparse_matrix();

  t_base       = uhm::timer();
  petsc->finalize();
  t_finalize = uhm::timer() - t_base;

  printf("==== Report =====\n");
  printf("Number of RHS          = %d\n", n_rhs);
  printf("Number of threads      = %d\n", uhm::get_num_threads());
  printf("Time init (s)          = %E\n", t_init);
  printf("Time solve (s)         = %E\n", t_solve);
  printf("Time check (s)         = %E\n", t_check);
  printf("Time finalize (s)      = %E\n", t_finalize);
  printf("--------------------------\n");


  delete petsc;

  m->unlock();

  delete m;

  FLA_Finalize();
  return 0;
}


